#!/usr/bin/perl -w
use strict;
use DBI;
use DBD::Pg;
use FindBin;
use lib ("$FindBin::Bin/..", "/net/cpp-group/Leo/bin");
use chromosomes;
use linewrap;
use uniquify;
use GO_frequency;
use Data::Dump qw(dump);

use db_parameters;
use ortho_species;

print STDERR <<"HEADLINE";
888888888888888888888888888888888888888888888888888888888888888888888

    table_significant_GO_terms

888888888888888888888888888888888888888888888888888888888888888888888
HEADLINE

#_________________________________________________________________________________________

#	get_duplicated_genes

#_________________________________________________________________________________________
sub get_all_gene_ids($$)
{
	my ($dbh, $species) = @_;
	print STDERR ("\tRetrieving all genes for $species...\n");

	my $sql_cmd = <<"PL/SQLCMD";
    SELECT
            gene_id
        FROM
            orthologs.predicted_type
        WHERE
            species = '$species' and
            protocol_id = $curr_protocol_id;
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd);
	print "\t\t", scalar @$data, " genes\n";
	#return uniquify(map{$_->[0] =~ /(ENS[A-Z]*G\d+)/; $1} @$data);
	return map{$_->[0]} @$data;
}


#_________________________________________________________________________________________

#	get_duplicated_genes

#_________________________________________________________________________________________
sub get_duplicate_genes($$$)
{
	my ($dbh, $species, $protocol_id) = @_;
	print STDERR ("\t\tRetrieving duplicates for $species...\n");

	my $sql_cmd = <<"PL/SQLCMD";
    SELECT
		  DISTINCT ON (ortholog_id, gene_id)
		  ortholog_id,
		  gene_id
	FROM
        orthologs.ortholog_sets
    WHERE
        protocol_id = $protocol_id AND
        ortholog_type ~ 'to' AND
        species = '$species'
	ORDER BY ortholog_id;
PL/SQLCMD
	my %all_clades;
	for my $data(@{$dbh->selectall_arrayref($sql_cmd)})
	{
		push(@{$all_clades{$data->[0]}}, $data->[1]);
	}
	my @genes;
	for my $ortholog_id(keys %all_clades)
	{
		# ignore clades with just one gene
		next unless @{$all_clades{$ortholog_id}} > 1;
		push(@genes, @{$all_clades{$ortholog_id}});
	}
	return uniquify(map{s/&.*//;$_} @genes);
	
}

#_________________________________________________________________________________________

#	get_pseudo_genes

#_________________________________________________________________________________________
sub get_pseudo_genes($$$)
{
	my ($dbh, $species, $protocol_id) = @_;
	print STDERR ("\t\tRetrieving pseudo genes for $species...\n");

	my $sql_cmd = <<"PL/SQLCMD";
    SELECT
		  gene_id
	FROM
        orthologs.ortholog_sets
    WHERE
        protocol_id = $protocol_id AND
        ortholog_type = 'pseudo' AND
        species = '$species';
PL/SQLCMD

	my $data = $dbh->selectall_arrayref($sql_cmd);
	
	my @data = uniquify(map{$_ = $_->[0]; s/&.*//; $_} @$data);
	print "\t\t", scalar @data, " genes\n";
	return @data;
}

#_________________________________________________________________________________________

#	get_homo_orphans

#_________________________________________________________________________________________
sub get_orphan_genes($$$)
{
	my ($dbh, $species, $protocol_id) = @_;
	print STDERR ("\t\tRetrieving $species orphaned gene ids...\n");

	my $sql_cmd = <<"PL/SQLCMD";
    SELECT
		DISTINCT ON (gene_id)
		gene_id
	FROM
		orthologs.orphans
	WHERE
		protocol_id = $protocol_id AND
        species = '$species';
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd);
	my @data = uniquify(map{$_ = $_->[0]; s/&.*//;$_} @$data);
}





#_________________________________________________________________________________________
#
#
#	retrieve_GO_ids_per_category
#
#_________________________________________________________________________________________
sub retrieve_GO_ids_per_category($\%$)
{
	print STDERR "\t\tGet GO identifiers per category per gene_id...\n";
	my ($dbh, $gene_id_to_GO_ids, $species) = @_;
	my $sqlcmd = <<"PL/SQLCMD";
    SELECT
        GO_category,
		gene_id,
		GO_id
	FROM 
		taxon.gene_go_matches JOIN
        GO_term_descriptions USING (GO_id)
    WHERE
        species = '$species';
PL/SQLCMD
	my $sql_data = $dbh->selectall_arrayref($sqlcmd);

	my $cnt = 0;
	# save GO identifiers per GO category per gene id
	for my $line(@$sql_data)
	{
		++$cnt;
		push(@{$gene_id_to_GO_ids->{$line->[0]}{$line->[1]}}, $line->[2]);
	}
	print "\t\t", $cnt, " GO terms\n";

}


#_________________________________________________________________________________________
#
#
#	retrieve_GO_descriptions
#
#_________________________________________________________________________________________
sub retrieve_GO_descriptions($\%)
{
	print STDERR "\tGet GO description...\n";
	my ($dbh, $GO_descriptions) = @_;
	my $sql_data = $dbh->selectall_arrayref(<<'PL/SQLCMD');
    SELECT
		GO_id, 
        description
	FROM 
		go_term_descriptions
PL/SQLCMD

	# save GO identifiers per GO category per gene id
	for my $line(@$sql_data)
	{
		$GO_descriptions->{$line->[0]} = $line->[1];
	}

}



#_________________________________________________________________________________________
#
#
#	retrieve_GO_ids_per_category
#
#_________________________________________________________________________________________
#sub retrieve_overrepresented_GO_ids($\%)
#{
#	print STDERR "\tGet significantly overrepresented GO IDs...\n";
#	my ($dbh, $overrepresented_GO) = @_;
#	my $sql_data = $dbh->selectall_arrayref(<<'PL/SQLCMD');
#    SELECT
#        protocol_id,
#		GO_id
#	FROM 
#		orthologs.significant_GO_terms
#    WHERE
#        category = 'duplicate'
#PL/SQLCMD
#
#	# save GO identifiers per GO category per gene id
#	for my $line(@$sql_data)
#	{
#		++$overrepresented_GO->{$line->[0]}{$line->[1]};
#	}
#
#}





#_________________________________________________________________________________________
#
#
#	print GO terms which are under or over represented
#
#_________________________________________________________________________________________
sub print_GO_over_under($$$$$\%\%\%\@\@)
{
	my ($dbh, $category, $protocol_id, $protocol_name, $species,
		$GO_slim_ids,
		$GO_id_to_description, $gene_id_to_GO_ids_all,
		$gene_ids, $gene_ids_background) = @_;

	print STDERR "\t\t\tGet GO over/underrepresentation...\n";

	my @GO_terms_type = (["go_slims", $GO_slim_ids, 1, 1e-2],
						 ["full_go", undef, 10, 1e-4]);

	for my $GO_terms_type(@GO_terms_type)
	{
		my ($GO_terms_type_name, $GO_id_filter, 
			$minimum_background_match_count, $p_value_threshold) = @$GO_terms_type;
		my $file_name = "$protocol_name.$species.$category.$GO_terms_type_name.go_cat";
		print STDERR "\t\t\t$GO_terms_type_name\n";

		# get significant results
		my (%significant_results, %parameters);
		get_under_overrepresented_GO_terms(	%$gene_id_to_GO_ids_all, $GO_id_filter, 
											 @$gene_ids_background, @$gene_ids, 
											 $p_value_threshold, 
											 %significant_results, %parameters);
	
	
		print STDERR "\t\t\tGO terms +/- calculated.\n";
		open (FILE, ">$dir_orthologs_output/go_terms/$file_name") or die;
		print FILE  scalar @$gene_ids, 				"\tgenes in sample population\n";
		print FILE scalar @$gene_ids_background,	"\tgenes in total population \n";
		for my $GO_cat (sort keys %significant_results)
		{
			print FILE "$parameters{$GO_cat}{sample_gene_id_pop}\tgenes in sample population with $GO_cat go terms.\n";
			print FILE "$parameters{$GO_cat}{tot_gene_id_pop}\tgenes in total population  with $GO_cat go terms.\n";
			print FILE "$parameters{$GO_cat}{tot_GO_id}\t$GO_cat GO terms in total population.\n";
		}
	
		
		# for each set of under and then over represented
		for my $over_under (qw(+ -))
		{
			# print out for each GO category
			for my $GO_cat (sort keys %significant_results)
			{
				# sort GO identifiers by their descriptions
				my @GO_ids = sort {	$GO_id_to_description->{$a} cmp
									$GO_id_to_description->{$b}}
									keys %{$significant_results{$GO_cat}{$over_under}};
	
				# for each GO identifier
				for my $GO_id(@GO_ids)
				{
					next unless $significant_results{$GO_cat}{$over_under}{$GO_id}[2] >
										$minimum_background_match_count;
					
					printf FILE	"$GO_cat\t$over_under\t$GO_id\t$GO_id_to_description->{$GO_id}\t".
									"%d\t%d\t%d\t%d\t%1.2g\t%1.2g\n",
									@{$significant_results{$GO_cat}{$over_under}{$GO_id}};
				}
			}
			
		}
	
		my %over_under_type = (	'+' => ['true', 4], #likelihoodG'],
								'-' => ['false', 5]); #'likelihoodF']		);
		
		print STDERR "\t\t\tCopying over/underrepresented GO terms IDs to panda\n";
		my $use_GO_slims = $GO_id_filter ? 1 : 0;
		$dbh->do("COPY orthologs.significant_GO_terms(over_representation, category, ".
				 "go_id, protocol_id, significance, species, GO_slim) FROM STDIN");
		# for each set of under and then over represented
		for my $over_under (qw(+ -))
		{
			my $over_bool = $over_under_type{$over_under}[0];
			my $prob_type = $over_under_type{$over_under}[1];
			# print out for each GO category
			for my $GO_cat (keys %significant_results)
			{
				my $results =  $significant_results{$GO_cat}{$over_under};
				# for each GO identifier
				for my $GO_id(keys %{$significant_results{$GO_cat}{$over_under}})
				{
					my $significance = $results->{$GO_id}[$prob_type];
					$significance = 0 if ($significance < 1e-300);
					$dbh->pg_putline("$over_bool\t$category\t$GO_id\t".
									"$protocol_id\t".$significance.
									"\t$species\t$use_GO_slims\n");
				}
			}
	
		}
		$dbh->pg_putline("\\.\n");
		$dbh->pg_endcopy;
		print STDERR "\t\t\tDone\n";
	}
}

#
#	Connect to genome database
#
my $dbh = connect_to_panda();
use get_psql_time;
my $start_time = get_psql_time($dbh);
#my @db_param = get_db_parameters();


#
#	Make output directory
#
system ("mkdir -p $dir_orthologs_output/go_terms");

#
#	GO identifiers and descriptions
#
my %GO_id_to_description;
retrieve_GO_descriptions($dbh, %GO_id_to_description);

#
#	GO slim identifiers
#
print STDERR "\tGet GO slim identifiers\n";
my %GO_slim_ids;
get_GO_slims(%GO_slim_ids, "$FindBin::Bin/goslim_generic.obo");
	
#
#	Overrepresented GO identifiers for each protocol
#
#print STDERR "\tGet overrepresented GO identifiers\n";
#my %overrepresented_GOids
#retrieve_overrepresented_GO_ids($dbh, %overrepresented_GOids)

#
#	for each protocol of interest
#


$dbh->do("TRUNCATE orthologs.significant_GO_terms");
PER_SPECIES:
for my $species($ortho_name1, $ortho_name2)
{
	#
	#	All genes for background
	#
	my @gene_ids = get_all_gene_ids($dbh, $species);


	#
	#	All go terms for background
	#
	my %gene_id_to_GO_ids_all;
	retrieve_GO_ids_per_category($dbh, %gene_id_to_GO_ids_all, $species);

	
	for my $protocol (get_orthologs_protocols())
	{
		my ($protocol_name,  $protocol_id) = @$protocol;
		next unless $protocol_id == $curr_protocol_id;
	
		print STDERR "\n\t$species / $protocol_name:\n";
	
		
		#
		#	Duplicates genes
		#
		my @duplicate_genes = get_duplicate_genes($dbh, $species, $protocol_id);
		print_GO_over_under($dbh, "duplicate", 
							$protocol_id, $protocol_name, $species,
							%GO_slim_ids,
							%GO_id_to_description, %gene_id_to_GO_ids_all,
							@duplicate_genes, @gene_ids);
#		last PER_SPECIES; # DEBUG
	
		#
		#	Pseudo genes
		#
		my @pseudo_genes = get_pseudo_genes($dbh, $species, $protocol_id);
		print_GO_over_under($dbh, "pseudo", 
							$protocol_id, $protocol_name, $species,
							%GO_slim_ids,
							%GO_id_to_description, %gene_id_to_GO_ids_all,
							@pseudo_genes, @gene_ids);
		#
		#	orphaned genes
		#
		my @orphan_genes = get_orphan_genes($dbh, $species, $protocol_id);
		print_GO_over_under($dbh, "orphan", 
							$protocol_id, $protocol_name, $species,
							%GO_slim_ids,
							%GO_id_to_description, %gene_id_to_GO_ids_all,
							@orphan_genes, @gene_ids);
	}
}









# write out sequences with mulitple sizes

log_pipeline_stage($dbh, $curr_protocol_id, 43, 'table significant_go_terms', $start_time);
$dbh->disconnect();

